数学建模

 

编程知识点

模型库-scipy-拟合-线性最小二乘拟合(leastsq)

 

Wang Haihua

🚅 🚋😜 🚑 🚔

1. 数据拟合

在科学研究和工程应用中经常通过测量、采样、实验等方法获得各种数据。对一组已知数据点集，通过调整拟合函数（曲线）的参数，使该函数与已知数据点集相吻合，这个过程称为数据拟合，又称曲线拟合。

插值和拟合都是根据一组已知数据点，求变化规律和特征相似的近似曲线的过程。但是插值要求近似曲线完全经过所有的给定数据点，而拟合只要求近似曲线在整体上尽可能接近数据点，并反映数据的变化规律和发展趋势。因此插值可以看作是一种特殊的拟合，是要求误差函数为 0 的拟合。

2. 线性最小二乘拟合

2.1 线性最小二乘拟合函数说明

线性最小二乘拟合是最简单和最常用的拟合方法。scipy.optimize 工具箱中的 leastsq()、lsq_linear()，scipy.stats 工具箱中的 linregress()，都可以实现线性最小二乘拟合。

2.1.1 scipy.optimize.leastsq 函数说明

leastsq() 根据观测数据进行最小二乘拟合计算，只需要观测值与拟合函数值的误差函数和待定参数的初值，返回拟合函数中的待定参数，但不能提供参数估计的统计信息。

leastsq() 可以进行单变量或多变量线性最小二乘拟合，对变量进行预处理后也可以进行多项式函数拟合。

scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

主要参数：

• func：可调用的函数，描述拟合函数的函数值与观测值的误差，形式为 error(p,x,y)，具有一个或多个待定参数 p。误差函数的参数必须按照 (p,x,y) 的顺序排列，不能改变。

• x0：一维数组，待定参数的初值。

• args：元组，线性拟合时提供观测数据值 (xdata, ydata)，观测数据 xdata 可以是一维数组（单变量问题），也可以是多维数组（多变量问题）。

返回值：

• x：一维数组，待定参数的最小二乘估计值。

2.2 Python 例程：单变量线性拟合

程序说明：

• scipy.optimize.leastsq() 与 scipy.stats.linregress() 都可以进行单变量线性拟合。leastsq() 既可以用于单变量也可以用于多变量问题；linregress() 只能用于单变量问题，但可以给出很多参数估计的统计结果。

• leastsq() 要以子函数来定义观测值与拟合函数值的误差函数，例程中分别定义了拟合函数 fitfunc1(p, x) 与误差函数error1(p, x, y) ，是为了方便调用拟合函数计算拟合曲线在数据点的函数值。注意 p 为数组 。

• leastsq() 中误差函数的函数名可以任意定义，但误差函数的参数必须按照 (p,x,y) 的顺序排列，不能改变次序。

• leastsq() 中观测数据 (x, yObs) 是以动态参数 args 的方式进行传递的。这种处理方式非常独特，没有为什么， leastsq() 就是这样定义的。
• linregress() 只要将观测数据 (x,yObs) 作为参数，默认单变量线性拟合，不需要定义子函数。

leastsq() 与 linregress() 进行线性拟合，得到的参数估计结果是相同的。

2.3 案例1（单变量拟合）

# 1. 单变量线性拟合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.optimize import leastsq  # 导入 scipy 中的最小二乘法拟合工具
from scipy.stats import linregress  # 导入 scipy 中的线性回归工具

def fitfunc1(p, x):  # 定义拟合函数为直线
    p0, p1 = p  # 拟合函数的参数
    y = p0 + p1*x  # 拟合函数的表达式
    return y

def error1(p, x, y):  # 定义观测值与拟合函数值的误差函数
    err = fitfunc1(p,x) - y  # 误差
    return err

# 创建给定数据点集 (x,yObs)
p = [2.5, 1.5]  # y = p[0] + p[1] * x
x = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
y = p[0] + p[1] * x  # 理论值 y
np.random.seed(1)
yObs = y + np.random.randn(x.shape[-1])  # 生成带有噪声的观测数据
# print(x.shape, y.shape, yObs.shape)

# 由给定数据点集 (x,y) 求拟合函数的参数 pFit
p0 = [1, 1]  # 设置拟合函数的参数初值
pFit, info = leastsq(error1, p0, args=(x,yObs))  # 最小二乘法求拟合参数
print("Data fitting with Scipy.optimize.leastsq")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}\np[1] = {:.4f}".format(pFit[0], pFit[1]))

# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fitfunc1(pFit,x)

Data fitting with Scipy.optimize.leastsq
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528

# 比较：线性回归，可以返回斜率，截距，r 值，p 值，标准误差
slope, intercept, r_value, p_value, std = linregress(x, yObs)
print("\nLinear regress with Scipy.stats.linregress")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}".format(intercept))  # 输出截距 intercept
print("p[1] = {:.4f}".format(slope))  # 输出斜率 slope
print("r^2_value: {:.4f}".format(r_value**2))  # 输出 r^2 值
print("p_value: {:.4f}".format(p_value))  # 输出 p 值
print("std: {:.4f}".format(std))  # 输出标准差 std

Linear regress with Scipy.stats.linregress
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
r^2_value: 0.9521
p_value: 0.0000
std: 0.1161

# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.text(8,3,"youcans-xupt",color='gainsboro')
ax.set_title("Data fitting with linear least squares")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="fitting curve")
plt.legend(loc="best")

plt.show()

Data fitting with Scipy.optimize.leastsq
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528

Linear regress with Scipy.stats.linregress
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
r^2_value: 0.9521
p_value: 0.0000
std: 0.1161

2.4 案例2（多变量线性拟合）

程序说明：

• scipy.optimize.leastsq() 既可以用于单变量也可以用于多变量问题，本例程求解一个二元线性拟合问题：$$y = p[0] + p[1] * x1 + p[2] * x2$$

• leastsq()求解多变量问题的方法与单变量问题类似，以子函数 error2(p, x1, x2, y) 来定义观测值与拟合函数值的误差函数，以动态参数 args 的方式传递观测数据 (x, yObs) 。

• leastsq() 中误差函数的函数名可以任意定义，但误差函数的参数必须按照 (p,x1,x2,y) 的顺序排列，不能改变次序。
• scipy 只能做一元线性回归，例程中通过调用 statsmodels.api 进行多元线性回归，可以得到各种统计参数，供读者参考。

# 2. 多变量线性拟合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.optimize import leastsq  # 导入 scipy 中的最小二乘法工具

def fitfunc2(p, x1, x2):  # 定义拟合函数为直线
    p0, p1, p2 = p  # 拟合函数的参数
    y = p0 + p1*x1 + p2*x2  # 拟合函数的表达式
    return y

def error2(p, x1, x2, y):  # 定义观测值与拟合函数值的误差函数
    err = fitfunc2(p, x1, x2) - y  # 计算残差
    return err

# 创建给定数据点集 (x,yObs)
np.random.seed(1)
p = [2.5, 1.5, -0.5]  # y = p[0] + p[1] * x1 + p[2] * x2
x1 = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
x2 = np.array([0., 1.0, 1.5, 2.2, 4.8, 5.0, 6.3, 6.8, 7.1, 7.5, 8.0])
z = p[0] + p[1]*x1 + p[2]*x2  # 理论值 z
zObs = z + np.random.randn(x1.shape[-1])  # 生成带有噪声的观测数据
print(x1.shape, z.shape, zObs.shape)

(11,) (11,) (11,)

# 由给定数据点集 (x,z) 求拟合函数的参数 pFit
p0 = [1, 1, 1]  # 设置拟合函数的参数初值
pFit, info = leastsq(error2, p0, args=(x1,x2,zObs))  # 最小二乘法求拟合参数
print("Data fitting with Scipy.optimize.leastsq:")
print("z = p[0] + p[1]*x1 + p[1]*x2")
print("p[0]={:.4f}\np[1]={:.4f}\np[2]={:.4f}".format(pFit[0], pFit[1], pFit[2]))

# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
zFit = fitfunc2(pFit, x1, x2)

Data fitting with Scipy.optimize.leastsq:
z = p[0] + p[1]*x1 + p[1]*x2
p[0]=2.6463
p[1]=2.2410
p[2]=-1.3710

# 多元线性回归：最小二乘法(OLS)
import statsmodels.api as sm
x0 = np.ones(x1.shape[-1])  # 截距列 x0=[1,...1]
X = np.column_stack((x0, x1, x2))  # (nSample,3): [x0, x1, x2]
model = sm.OLS(zObs, X)  # 建立 OLS 模型: y = b0*x0 + b1*x1 + b2*x2 + e
results = model.fit()  # 返回模型拟合结果
zFit = results.fittedvalues  # 模型拟合的 y值

print(results.summary())  # 输出回归分析的摘要
print("\nOLS model: y = b0*x0 + b1*x1 + b2*x2")
print('Parameters: ', results.params)  # 输出：拟合模型的系数

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.922
Model:                            OLS   Adj. R-squared:                  0.902
Method:                 Least Squares   F-statistic:                     47.21
Date:                Fri, 08 Apr 2022   Prob (F-statistic):           3.72e-05
Time:                        21:45:55   Log-Likelihood:                -17.382
No. Observations:                  11   AIC:                             40.76
Df Residuals:                       8   BIC:                             41.96
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.6463      0.942      2.808      0.023       0.473       4.820
x1             2.2410      1.043      2.148      0.064      -0.165       4.647
x2            -1.3710      1.312     -1.045      0.326      -4.396       1.654
==============================================================================
Omnibus:                        0.439   Durbin-Watson:                   3.144
Prob(Omnibus):                  0.803   Jarque-Bera (JB):                0.488
Skew:                          -0.082   Prob(JB):                        0.783
Kurtosis:                       1.981   Cond. No.                         35.4
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

OLS model: y = b0*x0 + b1*x1 + b2*x2
Parameters:  [ 2.64628055  2.24100973 -1.37104475]

D:\software_install\anaconda\lib\site-packages\scipy\stats\stats.py:1541: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=11
  warnings.warn("kurtosistest only valid for n>=20 ... continuing "

参考文献

• youcans Python小白的数学建模课 https://www.zhihu.com/column/c_1381900867424075776